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Abstract. 

We investigate the distribution of zeros of the partition function of the two- and 
three-dimensional symmetric ±J Ising spin glasses on the complex field plane. We 
use the method to analytically implement the idea of numerical transfer matrix which 
provides us with the exact expression of the partition function as a polynomial of 
fugacity. The results show that zeros are distributed in a wide region in the complex 
field plane. Nevertheless we observe that zeros on the imaginary axis play dominant 
\ roles in the critical behaviour since zeros on the imaginary axis are in closer proximity 

to the real axis. We estimate the density of zeros on the imaginary axis by an 
(S| . importance-sampling Monte Carlo algorithm, which enables us to sample very rare 

events. Our result suggests that the density has an essential singularity at the origin. 
This observation is consistent with the existence of Griffiths singularities in the present 
systems. This is the first evidence for Griffiths singularities in spin glass systems in 
equilibrium. 

<N 

1. Introduction 

In order to shed a unique light on the problem of finite-dimensional spin glasses, we 
discuss the distribution of zeros of the partition function of the ± J Ising spin glass in 
the complex field plane, the Lee- Yang zeros pp. All Lee- Yang zeros of a ferromagnetic 
system are proved to lie on the imaginary axis in the complex field plane from the 
circle theorem [I]. However, locations of Lee- Yang zeros for spin glass models in 
general do not obey simple rules. It is therefore interesting and important to study the 
distribution numerically in order to extract useful information on phase transitions in 
finite-dimensional spin glasses from a point of view different from conventional numerical 
methods like Monte Carlo simulations. The first study along this line was due to Ozeki 
and Nishimori [2]. From data for relatively small systems, they found the tendency which 
may be consistent with the existence of the Griffiths singularity [3j. Our present work 
is to expand their study to much larger systems and carry out quantitative evaluations 
of the density of zeros. 
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In the diluted ferromagnet, the Griffiths singularity means that the free energy is a 
non-analytic function of the field at the origin in the temperature range T C <T < T^ nre \ 
where T c is the critical temperature and T c ( pure ) is that of the corresponding non-random 
system. This temperature range is called the Griffiths phase. This unusual phenomenon 
in the otherwise-paramagnetic phase is caused by rare events of bond configurations, 
large-size connected clusters generated with small but non- vanishing probabilities. This 
fact was connected with an essential singularity in the density of zeros at the origin [4] 
because large clusters tend to have zeros close to the origin at temperatures below Tj pure ^ . 
The Griffiths singularity in the diluted ferromagnet has also been studied through the 
behaviour of the inverse susceptibility analytically [5] and numerically [6] . 

On the spin glass problem, Randeria et al studied the dynamics and suggested 
the existence of the Griffiths phase also in the spin glass systems [7J. There are few 
studies of Lee- Yang zeros of spin glass models in equilibrium except the work by Ozeki 
and Nishimori [2] although there are some papers [HI [9] on the distribution of zeros in 
the complex temperature plane, Fisher zeros [10]. In the present paper we evaluate the 
density of Lee- Yang zeros of the spin glass systems by an importance-sampling Monte 
Carlo algorithm 0, [EH CE21 EH E] ■ Our results give evidence that the density of zeros 
has an essential singularity, which suggests the existence of the Griffiths singularity in 
spin glasses. 

The outline of this paper is as follows. In section 2 we explain the exact evaluation 
of the partition function by numerical transfer matrix. Then, we show the distribution 
of zeros of the ±J Ising model on the complex field plane in section 3. In section 4, 
we analyze zeros that have direct relevance to the phase transition and compare our 
investigation with known characteristics of phase diagrams. In section 5, we examine 
the density of zeros using the importance-sampling Monte Carlo algorithm which enables 
us to sample very rare events, and we discuss the existence of the Griffiths singularity 
in the ±J Ising model. We present our conclusion in section 6. 



2. Exact partition function 

We consider the ±J Ising model in an external magnetic field h on the square and 
the simple cubic lattices with cylindrical boundary conditions, free in one direction and 
periodic in the others, which are suited for the transfer matrix method because we take 
the trace of spins one by one in spiral order. The Hamiltonian is 

H = ~J2 J ij a i a j -h^Vii (!) 

(hO) i 

where denotes nearest neighbour pairs and Oi = ±1. The interaction is chosen 
as J (> 0) or — J with probability 1/2. Using the number of states Q {E, M) for given 
E = | J2(ij) (1 ~~ Jij&i&j) and M = \ J2i (1 + 0*)? the partition function is expressed in 
terms of a two- variable polynomial as 

Z(x,y) =Trexp(-/?W) (2) 
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N s N b 



= y 



M 



(3) 



M=0 E=0 



where y = exp (2/3 /i) (fugacity), x = exp (2/? J) , N s is the number of sites (assumed to 
be even), and AT& is the number of bonds on the lattice. 

We focus on the ^/-dependence of Z (x, y) for a fixed x, since we are interested only 
in the field term for the investigation of Lee- Yang zeros. Thus, we treat the partition 
function as a single- variable polynomial 



To locate Lee- Yang zeros for the present model, we have to evaluate the partition 
function for a given set of random interactions analytically and solve the equation 
Z(y) = for y. 

There are three methods to evaluate the partition function at least partially 
analytically using the idea of numerical transfer matrix [15l [[6] . 

The first one is to evaluate the partition function at N s + 1 different external fields 
yi = exp (2(3 hi) (i = 1, 2, • • • , N s + 1) by using the numerical transfer matrix method. 
Then, we solve the set of linear equations Em=o c m?A M = Z (yi) (i = 1, 2, • • • , N s + 1) 
for cm (M = 0, 1, • • • , N s ) [2j. In this way we know the values of all coefficients of the 
polynomial Z (y) = J2m c mV M and therefore we can solve Z(y) = 0. In this process, 
since both x and y are treated as numbers, we have to be careful to keep very high 
precisions to avoid rounding errors in solving Z(y) = 0. In addition, this method is 
not very efficient because it requires a repetition of essentially the same calculations. 
An advantage is that the required memory size is very small and is easily parallelized 
on the grid computer. We use this technique for systems larger than 16 x 16 (see the 
Appendix). 

The second method, called the canonical transfer matrix, is to numerically give 
analytic partition functions at a fixed temperature for all values of M p2] . This method 
evaluates cm (written as Q (M, y) in [T7j) for all M at once. The required memory size 
is N s + 1 times larger than in the first method, but this process is done with only one 
calculation. Thus, the second method is faster than the previous method on a single 
CPU, but attention needs to be paid to the precision since we treat x as a numerical 
value. 

The third method is called the microcanonical transfer matrix [T7j, which 
evaluates perfectly exact partition functions for all values of M and E using symbolic 
manipulations on the computer. This algorithm does not involve numerical errors. 
However, the memory requirement is times larger than in the second method. 

In consideration of the balance between the computational time and the required 
memory size, we mainly use the canonical transfer matrix method in the present work. 
This method accompanies rounding errors. The necessary precision depends on the ratio 
of the largest and the smallest coefficients, which increases very rapidly with the system 
size or with the temperature decrease. We have verified the reliability of this approach 
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Figure 1. Schematic diagram for the distribution of zeros on the complex field plane. 



for the ferromagnetic Ising model of sizes comparable to the spin glass system we are 
interested in as explained in the Appendix. 

The system size we calculated is up to 20 x 20 in two dimensions and 4x5x8 
in three dimensions, where the last digit is for the free boundary direction and the 
others are for the periodic boundary conditions. We investigated at least about 10, 000 
samples of random interactions for each system size. We solved the polynomial equation 
Z(y) = by mathematica™ by specifying the precision to appropriate values. 

3. Distribution of zeros 

The free energy for a system with quenched randomness is defined as the average over 
samples for all possible bond configurations. This implies that the partition function of 
a quenched system is regarded as being given by the product of partition functions of 
all possible configurations, 



where F denotes the quenched free energy and J denotes the set of random interactions. 
This fact allows us to plot zeros of all randomly-chosen samples simultaneously on the 
complex plane of logy = 2/3h. The temperature will be measured in units of J/ks- A 
schematic diagram for the distribution of zeros is presented in figure [U 

Figure [2] shows the temperature dependence of the distribution of zeros for systems 
of size 10 x 10 in two dimensions and 4 x 4 x 6 in three dimensions. Throughout this 
paper the range of the real axis is from —20 to 20 in two dimensions and from —30 to 
30 in three dimensions, and the imaginary axis ranges from to m. 

Generally, zeros of both dimensions approach the real axis with decreasing 
temperature. Nearest zeros to the origin lie on the imaginary axis for both dimensions 
and all temperatures. Let us call the location of such zeros the edge (see figured]). A clear 
difference between two and three dimensions appears at low temperature. Zeros lying 
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off the imaginary axis but close to the origin (to be called the body) form a wedged-like 
shape in two dimensions whereas the body looks relatively rounded in three dimensions. 

We next discuss the dependence of the distribution on system size. Figures [3] and 
@] show the system-size dependence at T = 0.5 from system size 4 x 4 to 16 x 16 in two 
dimensions (figure [3]) and from 3x3x2 to 4x4x6 in three dimensions (figured]). The 
width of the distribution in the transverse direction does not change significantly as the 
system size increases and is likely to stay almost the same for sufficiently large systems. 
Both the edge and the body move toward the real axis as the system size increases in 
both dimensions. Especially, the body (i.e. zeros off the imaginary axis) tends to move 
more significantly than the edge. In the thermodynamic limit, if the body reaches the 
real axis away from the origin, it may imply that an ordered phase (that may exist along 
the real axis at and around the origin) would survive until the field exceeds a finite value 
corresponding to the point where the zero (the body) reaches the real axis. 

To analyze the density more closely, we divide the area of —20 < Re (2(3h) < 20, 
< Im (2(3h) < 7i (two dimensions) or -30 < Re (2f3h) < 30, < Im (2(3h) < n (three 
dimensions) into 400 x 200 boxes and count the number of zeros in each box. Figures [5] 
and [6] show the resulting density plots in two and three dimensions, respectively. It is 
readily seen that the density is very high on (and around) the imaginary axis in all cases. 
Regions of low density are cleared to white as can be verified by comparison of figure 
[3] (left-bottom panel) and figure [5] (bottom panel), for example. Zeros in those regions 
are expected to have no influence on the system properties in the thermodynamic limit. 
A marked difference between two and three dimensions is that the former develops a 
very sharp wedge in the body of the density profile as the size increases whereas, in the 
latter, the body remains rounded at its bottom part as commented above already. In 
two dimensions, the edge approaches the origin while the body does not show such an 
approach to the real axis at a faster rate than the edge with size increase. Also in three 
dimensions, the body seems to approach the real axis, but the difference between two 
cases in figure [6] is not clear since the difference of size is small. 

4. Approach to the real axis 

Based on the qualitative observation in the previous section, we analyze our data 
quantitatively in this section. We assume that the behaviour of the nearest zero to 
the real axis of each sample determines the phase transition of the ±J model in the 
thermodynamic limit. We therefore calculate the average location of the nearest zero 
of each sample and analyze the system-size dependence of this average location. For 
two dimensions, we choose the temperatures as T = rj pure ^ = 2/log(v / 2 + 1) (— 2.269) 
and T = 0.5 (in the presumed Griffiths phase). In three dimensions, the investigated 
temperatures are T — T g — 1.1 (the spin glass transition temperature p[8]) and T = 0.5 
(in the spin glass phase). 

Figures [7| and M show the real part (right part of figures) and the imaginary part 
(left part of figures) of the average location of nearest zeros as functions of the inverse 
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Figure 2. Temperature dependence of the distribution of zeros on the complex 2/3 h 
plane. The system sizes are 10 x 10 (left figures) and 4x4x6 (right figures). S denotes 
the number of samples. 
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Figure 3. System-size dependence in two dimensions at T = 0.5. 
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4x4x6,5= 12083 



Figure 4. System-size dependence in three dimensions at T = 0.5. 
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Figure 5. Density plot of the distribution of zeros in two dimensions at T = 0.5. 
Colours represent the density: red (highest), green, light blue, purple and white 
(lowest) in this order. The system size is 6 x 6 (top) and 16 x 16 (bottom). 
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Figure 6. Density plot of the distribution of zeros in three dimensions at T = 0.5. 
The same colour code is used as in figure The system size is 3 x 3 x 3 (top) and 
4x4x6 (bottom). 
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of the linear size L(= N^ d ) in two and three dimensions 0. 

We find that the average location of nearest zeros on the imaginary axis is lower than 
those off the imaginary axis for both dimensions as has been seen in previous figures. 
We may safely conclude that, if there is a phase transition caused by an approach of 
zeros to the origin, zeros on the imaginary axis reach the origin before those off the 
imaginary axis. 

In two dimensions, the imaginary part of the average location seems to approach 
the origin linearly with the inverse of system size. Solid lines in figure [7| show the best 
fits to linear functions. The average location has clearly a non-vanishing imaginary part 
in the limit L —> oc for T = T^ pure \ The data for T = 0.5 is marginal: It is difficult 
to discern whether or not the average location reaches the origin in the thermodynamic 
limit. Analysis from a different standpoint will be presented in the next section. The 
real part of the average location of zeros off the imaginary axis seems to rapidly approach 
the origin as L —> oc. 

In three dimensions, we tried the fit of the points on the imaginary axis to L _M , 
because we expect the zeros to reach the origin in the thermodynamic limit due to the 
existence of the spin glass phase. Solid lines on the left parts in figure [8] represent the 
best fits to the power function, L" L33(1) for T = T g = 1.1 and L' 1 ^ for T = 0.5. 

The value /i = 1.33(1) for the critical point T — T g — 1.1 in three dimensions may 
not be readily identified with a critical exponent: In the pure system it is known that 
the edge approaches the origin as L~ Vh [19]. However, in the present random system, 
we evaluated the average location of the nearest zeros, for which we have no established 
scaling that relates /i = 1.33(1) with critical exponents. In two dimensions our system 
does not show conventional critical behaviour at T = T^ pure \ and the exponent \± — 1 
in figure [71 would not have direct relevance to critical exponents. 

It is not easy to draw a definite conclusion on the behaviour of zeros in the 
thermodynamic limit in three dimensions since the system size is small. Nevertheless 
the following story is not inconsistent with our data: At T = T g) if the average location 
of zeros reaches the real axis, it is likely to be only at the origin. At low temperatures, 
some people may wish to interpret the data that the system-size dependence of average 
zeros off the imaginary axis suggests the possible stability of the spin glass phase in 
the external field. If the imaginary part of the average of zeros off the imaginary axis 
approaches (triangles on the left bottom of figure [8]) and the real part approaches a 
finite value (triangles on the right bottom), the spin glass phase would be stable in the 
external field. However, the data on the right-bottom panel of figure [8] may instead be 
extrapolated to 0, rather than to a finite value, in the thermodynamic limit. If this is 
indeed the case, it might imply the absence of the AT line [20] in the three dimensional 
±J Ising model. We should anyway be very careful to draw a conclusion from the 
present data for small systems. 

1 /3 

\ See the Appendix for reasons to choose N s ' as the linear size in three dimensions. 
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Figure 7. Average locations of nearest zeros as functions of the inverse of the linear 
size in two dimensions for T = x^ pure ) (filled symbols) and T = 0.5 (open symbols). 
The averages are calculated for zeros on the imaginary axis, off the imaginary axis, and 
both of these taken into account. Solid lines represent the best fits to linear functions. 
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Figure 8. Average locations of nearest zeros as functions of the inverse of the linear 
size in three dimensions for T = T g = 1.1 (top figures) corresponding to the spin glass 
transition temperature and T = 0.5 (bottom figures). Solid lines represent best fits to 
power functions, L -133 at T = 1.1 and L -1,59 at T = 0.5. 
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Figure 9. Density of the nearest zero at T = 0.5 for the system size 8 x 8 in the 
two-dimensional ± J model. The open squares have been obtained by the importance 
sampling MC algorithm, and filled squares are by the simple sampling of 10 5 samples. 



5. Density of zeros on the imaginary axis 

In order to analyze the behaviour of zeros in a more detailed way, we study the density of 
zeros on the imaginary axis in two dimensions. We observe that zeros on the imaginary 
axis are in close proximity to the real axis in figures [2] to El In particular we are interested 
in whether or not the density shows an essential singularity at the origin as an indicator 
of the Griffiths singularity The density of zeros on the imaginary axis will be denoted 
as g (0), where 9 is defined by 2/3h = iO for the location of the edge, the nearest zero 
to the origin on the imaginary axis, of each sample. In order to compare the density of 
the ±J model with that of the diluted ferromagnet, we also estimated the density for 
the bond-diluted ferromagnet with the bond probability being 1/2 at T = 1.0, which is 
believed to show an essential singularity at the origin [4]. 

The density of nearest zeros falls very rapidly as shown in figure [9] both for the 
conventional simple sampling and for the importance sampling (to be explained now). 
Since we are interested in the behaviour of g (9) where its value is extremely small, 
a special care should be taken to reduce statistical errors in the estimate of g(9). In 
general, the evaluation of events that occur with exponentially small probabilities is very 
difficult if one employs a simple-sampling procedure. We therefore used the importance- 
sampling Monte Carlo algorithm as follows 0, [HI O [13l [14] . 

The basic idea is to generate bond configurations according to the Metropolis 
algorithm. The new bond configuration J f is generated from the current configuration 
J with the probability 



where 9 (J) is the location of the edge of a given bond configuration J and P (9) is a 
guiding function. The guiding function is ideally equal to g (9) because equation ([6]) 




(6) 
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then guarantees that the stationary distribution of the Markov chain with equation ([6]) 
becomes the inverse of the guiding function, 1/P(9). Therefore, configurations with 
smaller g (9) are generated with higher probabilities. The idea is that the new bond 
configuration with a small probability P (9 (J)) is encouraged to be chosen. In order to 
proceed without the prior knowledge of g (0), we first choose an appropriate function as 
P (9 (J)) and try an incremental improvement as follows. 

After an update, the nearest zero 9 {J') is calculated from the new bond 
configuration J 7 , and the histogram H (9) is incremented at this 9 as H (9) \— H (9) + 1. 
Then the histogram H (9) will reach the stationary distribution, 1/P(9) multiplied by 
the true distribution g (#), 

H{0)ocg{0)xJ-. (7) 

Thus, we can estimate g (9) as 

g(d)<xH(6)xP(d). (8) 

In our calculations, the new bond configuration J' is generated by flipping a single bond 
out of the current bond configuration J. The initial guiding function P (9) is generated 
from the simple-sampling algorithm for about 10 5 independent bond configurations. We 
repeat the above process to update P (9) := H (9) x P (9) until H (9) becomes nearly 
flat for the desired range, and finally we suppose that P (9) is close enough to g (9). The 
estimated autocorrelation time of our data is about a few Monte Carlo steps. However, 
we employ the data taken at every step because we observed no visible differences in 
the log scale by the change of steps to take data. 

We have obtained the data for g (9) of the ± J model and the diluted ferromagnet. 
Figure [10 shows the estimated density of the edge for L = 6, 8, 10, 12 in two dimensions at 
T = 0.5 for the ±J model and T = 1.0 for the diluted ferromagnet. These temperatures 
are below the transition points of the corresponding pure ferromagnets and therefore 
the systems lie in the Griffiths phase if any. With an essential singularity at the origin 
in mind, we analyze the data by the following function in consideration of the finite-size 
effect 

A(L) 



g (0, L) « exp 



(e-0 o (L)) A 



(9) 



where 9 G (L) is an adjustable parameter expected to vanish in thermodynamic limit % 
The best fits to equation ([9]) have been obtained by choosing A = 2.2(1) for the ±J 
model and A = 1.0(1) for the diluted ferromagnet as shown in figure For the 

infinite system, 9 (L) is observed to vanish and the density becomes 

g{9) ^exp[-A/9 A ] . (10) 

§ A fit to a power, g (6) ex (0 — 0o (L)) -A , gave ^-dependent A, an inconsistent result. 
|| Error estimates in these exponents are chosen to be on the safe side such that the error bars in figure 
ITT1 are well included by the margin of several times. 
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Figure 10. Density of the nearest zero for sizes 6x6 to 12 x 12 in two dimensions at 
T = 0.5 for the ±J model (left) and T = 1.0 for the diluted ferromagnet (right). 
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Figure 11. 
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go (L) exp 
and A = 



fits 

A 



-A(L)/(6-e (L)) 



to the exponential function g (0, L) = 
with A = 2.2(1) for the ±J model (left) 



1.0(1) for the diluted ferromagnet (right), where go (L) denotes the 



parameter for normalization. 



This function has an essential singularity at 9 = 0. If A = 1, it is the same function as 
suggested by Bray and Huifang for the diluted ferromagnet [4J. It is remarkable that 
the analytical result of [4j has been confirmed numerically and, in addition, a similar 
result has been observed for the ±J model. This fact suggests that a common physics 
is likely to underlie the behaviour, which supports the existence of the Griffiths phase 
also in the ±J model although connected clusters are not trivially defined in the spin 
glass system. 
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6. Conclusion 

In this paper we have evaluated the zeros of the partition function of the symmetric ± J 
model in two and three dimensions by the transfer matrix method. The results revealed 
several outstanding features of the distribution of zeros in the complex field plane. 

On the qualitative aspects, we found that the distribution has a real part on the 
complex field plane. Nearest zeros to the origin lie on the imaginary axis. This edge 
moves toward the real axis as the system size increases in both dimensions. For the 
problem of the existence of the Griffiths singularity, we estimated the density of the 
nearest zero on the imaginary axis for each sample using the importance-sampling Monte 
Carlo algorithm. The idea of this algorithm is to enhance the probability of events with 
low probabilities. As a consequence, we could observe the density g (9) « exp(— A/9 A ) 
for both the ± J model and the diluted ferromagnet for very small 6 and extremely small 
g (9). This function has an essential singularity at 9 = 0: Thus this result for the density 
of the ±J model is compatible with the existence of the Griffiths singularity similarly 
to the diluted ferromagnet albeit through a slightly different mechanism as indicated 
by the difference between A = 2.2(1) for the ±J model and A = 1.0(1) for the diluted 
ferromagnet. This is the first evidence for a Griffiths singularity in spin glass systems in 
equilibrium. Further work is in progress to clarify the temperature dependence of the 
density and other properties using the method of [21 J for example, which may give us 
useful information on the phase transition. 
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Appendix A. The nearest zero of the pure system 

We can confirm the precision of our calculations from the computation of zeros of pure 
systems with uniform Jy = J > 0, for which all zeros are on the imaginary axis by the 
circle theorem. It is important to check the precision in these calculations, because we 
use the canonical transfer matrix which treats x as numerical values in the partition 
function Z (x,y). In tables [ATI and lA2l we list the location of the edge, nearest zero, as 
functions of the system size in two and three dimensions, respectively. These results for 
L > 17 have been obtained by the method of linear equations. 
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Table Al. Locations of the edge of the pure system in two dimensions for L = 3 to 
20 at T = T c (pure) (left) and T = 0.5 (right). 
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065449847087066 


3 x 


4 


X 


5 





138347107456545 





052359877648501 


4 x 


4 


X 


4 





133611326102681 





049087385315331 


4 x 


4 


X 


5 





108662540409281 





039269908236389 


4 x 


4 


X 


6 





092979491978941 





032724923521501 


4 x 


5 


X 


6 





078908669455350 





026179938817203 


4 x 


5 


X 


7 





069174127961650 





022439947553280 


4 x 


5 


X 


8 





062199985131278 





019634954106283 



Table A2. Locations of the edge of the pure system in three dimensions for 
N s = 3x3x2 to 4x5x8 at T = T c (pure) (left) and T = 0.5 (right). We used 
T c (pure) = 4.51 [22]. 
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Figure Al. Locations of the Figure A2. Locations of the edge 

edge as a function of the system as a function of system size N^ 1 ^ 3 . 

size The solid line is The solid lines are fitting results as 

the extrapolation by the BST = CN~ Vh/S . 
algorithm. 



To analyze the data, it is useful to remember that the angle 9 is scaled as [I 

f L~ Vh for T = T(P ure ) 

(L~ d for T < T c (pure) 

where y h = 15/8 (d = 2) and y h ~ 2.48 (d = 3) [22]. The location of the edge 9 is 
expected to follow the scaling at T = rj pure ^ jTTJ 

9 (L) = CL~ Vh{L) (l + dL-" + C 2 L~ 2u; + •••), (A.2) 

where yn{L) is defined as 

log \e (L + 1) /e (L)) 

* (£) = - log [(£ + !)/£] ' <A ' 3) 
The exponent is also expanded as 

(L) = (l + C[L~ UJ + C' 2 L~^ + •••). (A.4) 

We extrapolate (L) to L — > oc using the BST algorithm with a; = 1 [23]. The data 
in table ED yields ^ = 1.87476 (T = Tj pure )) and 2.00000 (T = 0.5) in two dimensions 
(figure ED), very close to equation (IA.1I) . In three dimensions, it is not trivial to define 
the linear size L, since the systems are not regular hexahedrons. Figure IA2I presents 
the analyses by the fitting function 9 = CN~ yh ' 3 . Then we obtained y^ = 2.49(3) 
(T = T c (P ure )) and 3.00000 (T = 0.5). These results agree well with equation (TO) . 

For further analysis we tried to fit the data at T = 0.5 to 9 = CL~ d to obtain 
C = 3.14163 (d = 2) and 3.14159 (d = 3). We thus conclude 

9 = 7rL~ d (T < T c ( pure) ) , (A.5) 

with negligibly small higher order corrections. These results give us justifications of our 
precision and of the use of N^ 3 as the linear size in three dimensions. 
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